setwd(choose.dir())
  
load("Doubs.RData")  


# Simple transformation of an environmental variable ==============

range(env$slo)
# Log-transformation of the slope variable (y = ln(x))
# Compare histograms and boxplots of the raw and transformed values
par(mfrow = c(2, 2))
hist(env$slo, 
     col = "bisque",
     right = FALSE, main="Histogram of slope")

hist(log(env$slo), 
     col = "light green", 
     right = FALSE, 
     main = "Histogram of ln(slope)"
)
boxplot(env$slo, 
        col = "bisque", 
        main = "Boxplot of slope", 
        ylab = "slope"
)
boxplot(log(env$slo), 
        col = "light green", 
        main = "Boxplot of ln(slope)",
        ylab = "log(slope)"
)

?log

log(0+1)

# Standardization of all environmental variables ==================

library(vegan)
?decostand
# compositionality bias (e.g. from high-throughput sequencing). see clr and rclr and Gloor 2017 paper!


# Center and scale = standardize the variables (z-scores)
env.z <- decostand(env, "standardize")


(env$dfs[1]-mean(env$dfs))/sd(env$dfs)


apply(env.z, 2, mean)# check if means = 0
apply(env.z, 2, sd)	# check if standard deviations = 1


# Same standardization using the scale() function (which returns 
# a matrix)
env.z2 <- as.data.frame(scale(env))
summary(env.z==env.z2)



# TASK 1: inspect the numerical environmental parameter of the mites dataset (mite.env)
#Q: which parameter are numerical?
#Q: which parameter requires transformation to improve normality?
#Q: which transformation improves normality?



###  ASSOCIATION MEASURES

# Remove empty site 8
spe <- spe[-8, ]
env <- env[-8, ]
spa <- spa[-8, ]


### Q-mode dissimilarity matrices
?vegdist
?dist

#visual inspection
spe

par(mfrow=c(1,1))
image(as.matrix(spe))

# modification of image() to include row and column axis labels
image.real <- function(mat) { 
  mat <- t(mat)[,nrow(mat):1]
  image(mat, axes = FALSE, col = hcl.colors(15, palette="viridis"))
  axis(1, at = seq(0, 1, length = nrow(mat)), labels = rownames(mat), las=2)
  axis(2, at = seq(0, 1, length = ncol(mat)), labels = colnames(mat))
  box() 
}

image.real(as.matrix(spe))

library(gplots)
heatmap.2(as.matrix(spe), Colv = FALSE, Rowv =FALSE, dendrogram="none", trace="none")


## Q-mode dissimilarity and distance measures for
## (semi-)quantitative data


# Percentage difference (aka Bray-Curtis) dissimilarity matrix
# on raw species data. ".db" means "distance Bray".
spe.db <- vegdist(spe)	# method = "bray" (default)
spe.db
image.real(as.matrix(spe.db))

# custom function coldiss (Francois Gillet)
coldiss(spe.db, nc = 15, diag = FALSE)


# Percentage difference (aka Bray-Curtis) dissimilarity matrix
# on log-transformed abundances (log(x+1))
spe.dbln <- vegdist(log1p(spe))
spe.dbln
image.real(as.matrix(spe.dbln))


# Chord distance matrix
spe.norm <- decostand(spe, "nor")
spe.dc <- dist(spe.norm)
image.real(as.matrix(spe.dc))


# Hellinger distance matrix
spe.hel <- decostand(spe, "hel")
spe.dh <- dist(spe.hel)
image.real(as.matrix(spe.dh))



## Q-mode dissimilarity measures
## for binary data

# Jaccard dissimilarity matrix using function vegdist()
spe.dj <- vegdist(spe, "jaccard", binary = TRUE)
image.real(as.matrix(spe.dj))

# Jaccard dissimilarity matrix using function dist()
spe.dj2 <- dist(spe, "binary")



# Sorensen dissimilarity matrix
spe.ds <- vegdist(spe, method = "bray", binary = TRUE)
image.real(as.matrix(spe.ds))



## Compare distance matrices from environmental, species and
## spatial data

# Remove the 'dfs' variable from the env dataset
env2 <- env[, -1]

# Euclidean distance matrix of the standardized env2 data frame
env.de <- dist(scale(env2))
image.real(as.matrix(env.de))

coldiss(env.de, nc=8, diag=TRUE)

# Euclidean distance matrix on spatial coordinates (2D)
plot(spa$X, spa$Y, type="b", col="blue", xlab="X", ylab="y")
spa.de <- dist(spa)
coldiss(spa.de, nc = 16, diag = TRUE)

# Euclidean distance matrix on distance from the source (1D)
#dfs.df <- as.data.frame(env$dfs, row.names = rownames(env))
#head(dfs.df)
#riv.de <- dist(dfs.df)
#coldiss(riv.de, nc = 16, diag = TRUE)



# Fictitious data for Gower index
# Random normal deviates with zero mean and unit standard deviation
?rnorm
var.g1 <- rnorm(30, 0, 1)
var.g1
mean(var.g1)
hist(var.g1)

# Random uniform deviates from 0 to 5
var.g2 <- runif(30, 0, 5)
?runif
var.g2
hist(var.g2)


# Factor with 3 levels (10 objects each)
var.g3 <- gl(3, 10, labels = c("A", "B", "C"))
var.g3
# Factor with 2 levels, orthogonal to var.g3
var.g4 <- gl(2, 5, 30, labels = c("D", "E"))

(dat2 <- data.frame(var.g1, var.g2, var.g3, var.g4))
head(dat2)
summary(dat2)

# Computation of a matrix of Gower dissimilarity using
# function daisy()

# Complete data matrix (4 variables)
dat2.S15 <- daisy(dat2, "gower")
range(dat2.S15)
coldiss(dat2.S15, diag = TRUE)


#TASK 2: use the mite.spe and plot the Bray-Curtis dissimilarity and Jaccard dissimilariy matrices
#Q: visually estimate how many (large) clusters of samples are there based on Bray-Curtis dissimilarity?
#Q: is there a difference between number of clusters of samples based on Bray-Curtis and Jaccard dissimilarities?
#Q: Based on the results to the question above - what can we learn about the distribution of species and their abundance in the mite dataset?




### R-mode 


## R-mode correlation matrices

# Pearson r linear correlation among environmental variables
env.pearson <- cor(env)	# default method = "pearson"
round(env.pearson, 2)
coldiss(env.pearson)
pairs(env)  # see "improved" version in script 1


# Kendall tau rank correlation among environmental variables
env.ken <- cor(env, method = "kendall")
env.ken
env.o <- order.single(env.ken)
env.o

?order.single


panel.hist <- function(x, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(usr[1:2], 0, 1.5) )
  h <- hist(x, plot = FALSE)
  breaks <- h$breaks; nB <- length(breaks)
  y <- h$counts; y <- y/max(y)
  rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor*r)
}

pairs(
  env[, env.o],
  lower.panel = panel.smooth,
  upper.panel = panel.cor,
  no.col = TRUE,
  method = "kendall",
  diag.panel = panel.hist,
  main = "Kendall Correlation Matrix", cex.cor=2
)